#from scipy import *
#It is easy to not use the quad
#from scipy.integrate import quad
from numpy import *
from scipy.special import jn

# wavenumber
k = pi
#dipole pozition
z0 = 0

# koefficients for Gauss quadrature formula
d = array([0.0950125098376374401853193,0.2816035507792589132304605,0.4580167776572273863424194,0.6178762444026437484466718,0.7554044083550030338951012,0.8656312023878317438804679,0.9445750230732325760779884,0.9894009349916499325961542])
d1 = d[::-1]
d = hstack((-d1,d))
D = array([0.1894506104550684962853967,0.1826034150449235888667637,0.1691565193950025381893121,0.1495959888165767320815017,0.1246289712555338720524763,0.0951585116824927848099251,0.0622535239386478928628438,0.0271524594117540948517806])
D1 = D[::-1]
D = hstack((D1,D))

# Gauss quadrature formula
def quad_Gauss(f):
    return sum(D*f(d))


# contour parametrization and its derivatives [alpha,betta]->R
alpha = 0.0
betta = 1.0
def ro1(t):    return sin(t)
def z1(t):    return cos(t)
def dro1(t):    return cos(t)
def dz1(t):    return -sin(t)
def ddro1(t):    return -sin(t)
def ddz1(t):    return -cos(t)
def dddro1(t):    return -cos(t)
def dddz1(t):    return sin(t)



#length of contour
def len(t):
    y = alpha*(1-t)/2 + betta*(1+t)/2
    return sqrt(dro1(y)**2+dz1(y)**2)*(betta-alpha)/2 

length = quad_Gauss(len)

#discretization order
n = int(3*length*k/pi) + 5.
print "n =",n

# the  roots of Chebyshev polynomials(Chebyshev points)
ind1 = arange(int(n-1))
ind2 = arange(int(n))
t0 = cos(pi*(ind1+1)/n)
t = cos(pi*(2*ind2+1)/(2*n))
ind = arange(2*n-1)
t0t = hstack((t0,t))

# values of contour parametrization 
# and its derivatives [-1,1]->R in the Chebyshev points
mult = (betta-alpha)/2
tau = alpha*(1-t0t)/2 + betta*(1+t0t)/2
ro = ro1(tau)
z = z1(tau)
dz = dz1(tau)*mult
dro = dro1(tau)*mult
ddz = ddz1(tau)*mult**2
ddro = ddro1(tau)*(mult**2)
dddro = dddro1(tau)*(mult**3)
dddz = dddz1(tau)*(mult**3)

#Lyame koefficient
h_tau = sqrt(dro**2+dz**2)

#distance between to point on the surface
def L(m,p,psi): return sqrt(ro[p]**2+ro[m]**2-2*ro[p]*ro[m]*cos(psi)+(z[p]-z[m])**2)

#FUCKING LIMITS

def f_Smooth(M,p,psi):
    LL = sqrt(2)*ro[p]*sqrt(1-cos(psi))
    Part11 = -(dro[p]**2)*cos(psi)-dz[p]**2
    Part12 = exp(-1j*k*LL)/(LL**2)*(-1/LL - 1j*k)*cos(M*psi)
    Part13 = (1/LL**3)*(1 - (M**2 * sin(psi)**2)/2)
    Part14 = (k**2/2 + M**2/(2*ro[p]**2)*(cos(psi) - 1))/LL
    Part21 = (ro[p]*dro[p]*(1 - cos(psi)))**2
    Part22 = exp(-1j*k*LL)/(LL**3)*(-k**2 + 3*1j*k/LL +3/LL**2)*cos(M*psi)
    Part23 = -3*(1 - (M**2 * sin(psi)**2)/2)/(LL**5)
    Part24 = -(1/LL**3)*(k**2/2 + M**2/(2*ro[p]**2)*(cos(psi)-1))
    Res = Part11*(Part12 + Part13 + Part14) + Part21*(Part22+Part23+Part24)
    return Res
    
def Smooth(M,p):
    #number of oscilation in the function f_Smooth    
    num_oscilation = (2*ro[p]*k)/pi + M
    #Gauss quadrature formula with 16 points approximate function with 8 oscilations with accuracy 10^-7    
    # num is the number of intervals. For each interval we use Gauss formula with 16 points     
    num = int(num_oscilation/6) + 1  
    Int = 0    
    for i in xrange(num):
        arg = pi*(2*i + 1)/(2*num) + d*pi/(2*num)
        Int = Int + sum(D*f_Smooth(M,p,arg))    
    return Int*pi/num

def Lim_1_L(p):
    Part1 = (-dro[p]**2-dz[p]**2)/(2*ro[p]**2)
    mult1 = ddro[p]**2/2 + dro[p]*dddro[p]/3 + ddz[p]**2/2 + dz[p]*dddz[p]/3
    mult2 = dro[p]**2 + dz[p]**2
    mult3 = (dro[p]*ddro[p] + dz[p]*ddz[p])**2
    Part2 = -(mult1*mult2-mult3)/mult2**2
    Part3 = -dro[p]/ro[p]**2*(dro[p]/ro[p]-(dro[p]*ddro[p]+dz[p]*ddz[p])/mult2)
    Part4 = (dz[p]**2 + 2*dro[p]**2)/(2*ro[p]**3)*log(8*ro[p]/sqrt(mult2))    
    Part5 = mult2/(2*ro[p])**3*(5-2*log(8*ro[p]/sqrt(mult2)))   
    Part6 = -(3*dz[p]**2 + 2*dro[p]**2 + 2*ro[p]*ddro[p])/(4*ro[p]**3)   
    return (Part1 + Part2)/ro[p] + Part3 + Part4 + Part5+ Part6

def Lim_L(p):
    mult = dro[p]**2 + dz[p]**2    
    return -2*dz[p]**2/ro[p] + 2*mult/ro[p]*(2 - log(8*ro[p]/sqrt(mult)))

def Lim_ddS(M,p): 
    return Smooth(M,p) + Lim_1_L(p) + (M**2/(2*ro[p]**2) - k**2/2)*Lim_L(p)

def f_A(M,p,psi):
    LL = sqrt(2)*ro[p]*sqrt(1-cos(psi))    
    return (exp(-1j*k*LL)*cos(M*psi)-1)/LL

def Smooth2(M,p):   
    num_oscilation = (2*ro[p]*k)/pi + M
    num = int(num_oscilation/6) + 1  
    Int = 0      
    for i in xrange(num):
        arg = pi*(2*i + 1)/(2*num) + d*pi/(2*num)
        Int = Int + sum(D*f_A(M,p,arg))    
    return Int*pi/num    

def Lim_S(M,p): 
    return Smooth2(M,p) + log(8*ro[p]/sqrt(dro[p]**2 + dz[p]**2))*2/ro[p]

def Lim_K1tau(M,p):
    Part = dro[p]**2*(Lim_S(M+1,p)+Lim_S(M-1,p))/2 + dz[p]**2*Lim_S(M,p)    
    return ro[p]**3*(Lim_ddS(M,p) - k**2*Part)

def Lim_K2phi(M,p):
    return ro[p]**3*(M**2/ro[p]*Lim_S(M,p) - k**2*ro[p]/2*(Lim_S(M-1,p) + Lim_S(M+1,p)))

#quadrature formulas for the hypersingular, singular integrals and integrals with the logarithm kernels


def A(l,j): 
    if l == j: return -n/2
    else: return (1-(-1)**(j+l))*(1-t0[j]**2)/(t0[l]-t0[j])**2/n

def Gamma1(l,j):
    if l == j: return 0
    else: return (1-(-1)**(j+l))*(1-t0[j]**2)/(t0[l]-t0[j])/n

def T(p,t): 
    return cos(p*arccos(t))

def LI(tt,j):
    part = (T(ind1+1,t0[j])*T(ind1+1,tt)/(ind1+1)).sum()
    return -1/n*(1-t0[j]**2)*(log(2) + 2*part + (-1)**j*T(n,tt)/n)

def LII(tt,j):
    part = (T(ind1+1,t[j])*T(ind1+1,tt)/(ind1+1)).sum()
    return -1/n*(log(2) + 2*part)

def Gamma2(l,j): 
    return (1/n)*1/(t0[l]-t[j])

def Gamma3(l,j): 
    return (1/n)*(1-t0[j]**2)/(t[l]-t0[j])

#quadrature formulas for smooth kernels
#simplify
def m0(M,p,j):
    Lmin = sqrt((ro[j]-ro[p])**2 + (z[j]-z[p])**2)
    Lmax = sqrt((ro[j]+ro[p])**2 + (z[j]-z[p])**2)
    a = (Lmax - Lmin)/(3*Lmin + Lmax)   
    arg = 10**(-4)*(1-a**2)/pi            
    nn = abs(log10(arg)/log10(a)) + 2
    return round((Lmax-Lmin)*k/pi + nn + M)


def a1_tau(p): 
    return -2*ro[p]**2

def b1_tau(p): 
    return -dro[p]*ro[p]

def c1_tau(M,p): 
    return (k**2*(ro[p])**2 + M**2)*(dro[p]**2 + dz[p]**2) - (dz[p]**2 + 3*dro[p]**2)/4

def dg(p,j,psi): 
    return ro[j]*dro[j] - ro[p]*dro[j]*cos(psi) + (z[j] - z[p])*dz[j]

def d2g(p,j,psi):
    return -dro[p]*dro[j]*cos(psi) - dz[p]*dz[j]

def f_A_1tau(M,p,j,psi):
    LL = L(p,j,psi)
    Part1 = d2g(p,j,psi)*(1/LL**2)*(-1/LL - 1j*k)
    Part2 = dg(p,j,psi)*dg(j,p,psi)*1/LL**3*(-k**2 + 3*1j*k/LL + 3/LL**2)
    Part3 = -k**2*(dro[p]*dro[j]*cos(psi)/LL + dz[p]*dz[j]/LL)
    return exp(-1j*k*LL)*(Part1 + Part2 + Part3)*cos(M*psi)


#PROBLEM - IT IS WRONG ARG ARANGE(m), if m IS VECTOR

# quadrature formula in which function defined in [0,pi] interpolate by trigonometric polynomials
def Quad_trig(f,M,p,j):
    m = int(m0(M,p,j))
    phi = 2*pi*(arange(m)+1)/(2*m+1)      
    part = sum(f(M,p,j,phi))   
    return 2*pi/(2*m+1)*(f(M,p,j,0) + 2*part)

def K1_tau(M,p,j): 
    if p!=j:    
        Part1 = ro[p]**3*Quad_trig(f_A_1tau,M,p,j) - a1_tau(p)/(t0[p]-t0[j])**2 
        Part2 = - b1_tau(p)/(t0[p]-t0[j]) - c1_tau(M,p)*log(abs(t0[p]-t0[j]))       
        return (Part1 + Part2)*(1-t0[j]**2)/n
    else:
        return Lim_K1tau(M,p)*(1-t0[p]**2)/n

def b1_phi(M,p):
    return 2*1j*M*ro[p]**2

def c1_phi(M,p):
    return -1j*M*dro[p]*ro[p]

def f_A_1phi(M,p,j,psi):
    LL = L(p,j,psi)
    Part1 = (-1j*M)*(ro[p]*dro[p] - dro[p]*ro[j]*cos(psi) + (z[p] - z[j])*dz[p])/LL**2*(-1/LL - 1j*k)
    Part2 = k**2*ro[j]*dro[p]*(cos((M+1)*psi) - cos((M-1)*psi))/(2*1j*LL)
    return exp(-1j*k*LL)*(Part1*cos(M*psi) + Part2)

def K1_phi(M,p,j): 
    Part1 = ro[p]**3*Quad_trig(f_A_1phi,M,p,j + n-1) - b1_phi(M,p)/(t0[p]-t[j])
    Part2 = -c1_phi(M,p)*log(abs(t0[p]-t[j]))       
    return (Part1 + Part2)/n

def b2_tau(M,p):
    return 2*1j*M*ro[p]

def c2_tau(M,p):
    return 1j*M*dro[p]

def f_A_2tau(M,p,j,psi):
    LL = L(p,j,psi)
    Part1 = (1j*M)/ro[p]*(ro[j]*dro[j] - dro[j]*ro[p]*cos(psi) + (z[j] - z[p])*dz[j])/LL**2*(-1/LL - 1j*k)
    Part2 = -k**2*dro[j]*(cos((M+1)*psi) - cos((M-1)*psi))/(2*1j*LL)
    return exp(-1j*k*LL)*(Part1*cos(M*psi) + Part2)

def K2_tau(M,p,j): 
    Part1 = ro[p+n-1]**3*Quad_trig(f_A_2tau,M,p+n-1,j) - b2_tau(M,p + n-1)/(t[p]-t0[j])
    Part2 = -c2_tau(M,p+n-1)*log(abs(t[p]-t0[j]))       
    return (Part1 + Part2)*(1-t0[j]**2)/n

def c2_phi(M,p):
    return -2*M**2*ro[p] + 2*k**2*ro[p]**3

def f_A_2phi(M,p,j,psi):
    LL = L(p,j,psi)
    Part = M**2/ro[p]/LL - k**2*ro[j]/LL*cos(psi)
    return exp(-1j*k*LL)*Part*cos(M*psi)

def K2_phi(M,p,j): 
    if p!=j:
        Part = -c2_phi(M,p+n-1)*log(abs(t[p]-t[j]))       
        return (ro[p+n-1]**3*(Quad_trig(f_A_2phi,M,p+n-1,j+n-1)) + Part)/n
    else:
        return Lim_K2phi(M,p+n-1)/n




#MAIN MATRIX

def Matr(M):
    C1_tau = zeros((n-1,n-1),complex)
    C1_phi = zeros((n-1,n),complex)
    C2_tau = zeros((n,n-1),complex)
    C2_phi = zeros((n,n),complex)
    for p in xrange(0,int(n)-1):
        for j in xrange(0,int(n)-1):
            C1_tau[p,j] = a1_tau(p)*A(p,j) + b1_tau(p)*Gamma1(p,j) + c1_tau(M,p)*LI(t0[p],j) + K1_tau(M,p,j)
    for p in xrange(0,int(n)-1):
        for j in xrange(0,int(n)):
            C1_phi[p,j] = b1_phi(M,p)*Gamma2(p,j) + c1_phi(M,p)*LII(t0[p],j) + K1_phi(M,p,j)
    for p in xrange(0,int(n)):
        for j in xrange(0,int(n)-1):
            C2_tau[p,j] = b2_tau(M,p+n-1)*Gamma3(p,j) + c2_tau(M,p+n-1)*LI(t[p],j) + K2_tau(M,p,j)                          
    for p in xrange(0,int(n)):
        for j in xrange(0,int(n)):
            C2_phi[p,j] = c2_phi(M,p+n-1)*LII(t[p],j) + K2_phi(M,p,j)
    return vstack((hstack((C1_tau,C1_phi)), hstack((C2_tau,C2_phi))))                                                            
                          
def R(r,z):
    return sqrt(r**2 + z**2)

# electrical dipole electric field components
def E0_R_0(r,z):
    R1 = R(r,z)
    part = -1j/k*(1/R1**2 + 1j*k/R1)
    return 30*exp(-1j*k*R1)/R1*(r/sqrt(r**2+z**2)*part)

def E0_teta_0(r,z):
    R1 = R(r,z)
    part = 1j/k*(1/R1**2 + 1j*k/R1 - k**2)
    return 30*exp(-1j*k*R1)/R1*(z/sqrt(r**2+z**2)*part)*0.5

def E0_phi_0(r,z):
    R1 = R(r,z)
    part = 1j/k*(1/R1**2 + 1j*k/R1 - k**2)
    return 30*exp(-1j*k*R1)/R1*part*0.5j

def E0_r_0(r,z):
    return r/sqrt(r**2+z**2)*E0_R_0(r,z) + z/sqrt(r**2+z**2)*E0_teta_0(r,z) 

def E0_z_0(r,z):
    return z/sqrt(r**2+z**2)*E0_R_0(r,z) - r/sqrt(r**2+z**2)*E0_teta_0(r,z) 

# the right side of the matrix equation

f_tau = 4*1j*k*(dro*E0_r_0(ro, z-z0) + dz*E0_z_0(ro, z-z0))/(120*pi)*ro**3
f_phi = 4*1j*k*E0_phi_0(ro, z-z0)/(120*pi)*ro**3

f_tau = f_tau[0:n-1]
f_phi = f_phi[n-1:2*n-1]


f = hstack((f_tau,f_phi))

#solution of the matrix equation
M = Matr(1)

x = linalg.solve(M, f)

x_tau = x[0:n-1]
x_phi = x[n-1:2*n-1]

def f(M,teta,p):
    return 1j**M*2*pi*exp(1j*k*cos(teta)*z[p])*jn(M,k*ro[p]*sin(teta))

#print f(1,1,n-1+1)

n = int(n)
def F_phi(M,teta):
    vect_tau = k**2*dro[0:n-1]*(f(M+1,teta,arange(n-1)) - f(M-1,teta,arange(n-1)))/2j*(1-t0[0:n-1]**2)*x_tau
    vect_phi = k**2*ro[n-1:2*n-1]*(f(M+1,teta,arange(n)+ n-1) + f(M-1,teta,arange(n) + n-1))/2*x_phi   
    return (-1j)*30*pi/k*(sum(vect_tau)+sum(vect_phi))/n


def F_teta(M,teta):
    vect_tau = (dz[0:n-1]*f(M,teta,arange(n-1))*sin(teta)-dro[0:n-1]*(f(M+1,teta,arange(n-1)) + f(M-1,teta,arange(n-1)))*cos(teta)/2)*(1-t0[0:n-1]**2)*x_tau
    vect_phi = ro[n-1:2*n-1]*(f(M+1,teta,arange(n)+ n-1) - f(M-1,teta,arange(n) + n-1))/2j*cos(teta)*x_phi   
    return 1j*30*pi*k*(sum(vect_tau)+sum(vect_phi))/n

print F_teta(1,1)


from time import time
t1 = time()
t2 = time()
print t2-t1





